#include "macros.hpp"
#include "eigenHelpers.hpp"

namespace hms::swe
{

/* this should be implemented in RectCells, because it needn't be recalculated */
template<typename MeshType, EdgeRange orien, int side>
HMS_INLINE auto dxQuotient(
	[[maybe_unused]] const Ref<const ArrayXs>& dx
){
	if constexpr ( std::is_same_v<MeshType, UniMesh> ){
		return 1;
	} else if constexpr ( std::is_same_v<MeshType, RectMesh>){
		if constexpr ( orien == EdgeRange::horizontal ){
			/* now dx is actually dy, and dxQuotient is scalar */
			assert( dx.size() == 3 );
			return dx[1] / dx[2*side];
		} else {
			/* only for vertical edges, dxQuotient is an array */
			Index nEdges { dx.size() - 2 };
			assert( dx.size() == nEdges + 2 );
			return (
				dx.segment(1, nEdges) /
				( side == 0 ? dx.head(nEdges) : dx.tail(nEdges) )
				)
				.reshaped().transpose(); // that always creates a row vector
		}
	}
}


template<typename SolverType>
void murilloCheck(
	EdgesValues& ev,
	const Pair<Array1XsCMap>& dCellLR,
	const Pair<Array1XsCMap>& zCellLR,
	const Pair<Array2XsCMap>& vCellLR
){
	using S = SolverType;

	Index nEdges { dCellLR[0].cols() };

	/* we currently use a loop for two reasons:
	 * 1. the murillo condition only has to be evaluated once for both sides,
	 * which didn't perform as well using replicate<1,2>()
	 * (likely because it branches twice as often)
	 * 2. firstOrder() is implemented as a void function, which means it would
	 * have to get evaluated for all values beforehand
	 * */
	scalar dL, dR;
	for (int i{0}; i<nEdges; ++i){
		dL = dCellLR[0][i];
		dR = dCellLR[1][i];
		/* this is not only Murillo's condition, but also the check whether
		 * one of the two water depths is below zero.
		 * In both cases, first order reconstruction is performed, so they're
		 * combined here.
		  */
		if (
			min( dL, dR ) < max( max(dL, dR) - min(dL, dR), S::minWaterDepth )
		){
			for (int side{0}; side<2; ++side){
				reconstruction::firstOrder<S>(
					ev.dLR[side][i], ev.zLR[side][i], ev.hLR[side][i], ev.vLR[side].col(i),
					dCellLR[side][i], zCellLR[side][i], vCellLR[side].col(i)
				);
			}
		}
	}
}


// template<typename SolverType, EdgeRange orien, typename A, typename B>
// void secondOrderReconstruction_qh(
// 	EdgesValues& ev,
// 	const Pair<A>& dCellLR,
// 	const Pair<A>& dCellLLRR,
// 	const Pair<A>& hCellLR,
// 	const Pair<A>& hCellLLRR,
// 	const Pair<B>& qCellLR,
// 	const Pair<B>& qCellLLRR,
// 	[[maybe_unused]] const Ref<const ArrayXs>& dx,
// 	const LimiterFunction& limFunc
// ){
// 	static_assert( std::is_base_of_v<SolverBase<SolverType>, SolverType> );
// 	using S = SolverType;
// 	using MeshType = typename S::MeshType;

// 	ASSERT_MATCHING_SIZES(dCellLR, hCellLR, qCellLR);

// 	assert( ev.dLR[1].data() == ev.dLR[0].data() + ev.dLR[0].size() );
// 	assert( ev.hLR[1].data() == ev.hLR[0].data() + ev.hLR[0].size() );
// 	assert( ev.vLR[1].data() == ev.vLR[0].data() + ev.vLR[0].size() );

// 	auto
// 		dxL { dxQuotient<MeshType, orien, 0>(dx) },
// 		dxR { dxQuotient<MeshType, orien, 1>(dx) };

// 	reconstruction::highOrder<S>(
// 		ev,
// 		dCellLR, dCellLLRR, hCellLR, hCellLLRR, qCellLR, qCellLLRR,
// 		dxL, dxR, limFunc
// 	);
// }





// template<typename T>
// struct MapHelper {};

// template<typename PlainObjectType, int MapOptions, typename StrideType>
// struct MapHelper<Map<PlainObjectType, MapOptions, StrideType>> {
// 	using Type = PlainObjectType;
// };

// template<typename T>
// using MapType = typename MapHelper<T>::Type;


template<typename MeshType, EdgeRange orien>
void transformVelocities( EdgesValues& ev ){

	#ifndef NDEBUG
	Index nEdges { ev.vLR[0].cols() };
	assert( nEdges == ev.vLR[1].cols() );
	assert( nEdges == ev.vRotLR[0].cols() );
	assert( nEdges == ev.vRotLR[1].cols() );
	#endif

	for (int i{0}; i<2; ++i){
		if constexpr (
			   std::is_same_v<MeshType, UniMesh>
			|| std::is_same_v<MeshType, RectMesh>
		){
			if constexpr ( orien == EdgeRange::horizontal ){
				ev.vRotLR[i].row(0) =  ev.vLR[i].row(1);
				ev.vRotLR[i].row(1) = -ev.vLR[i].row(0);
			} else {
				ev.vRotLR[i] = ev.vLR[i];
			}
		} else {
			assert(false);
		}
	}
}

} // namespace hms::swe

#undef ASSERT_MATCHING_SIZES